﻿/*
* Copyright 1993-2015 NVIDIA Corporation.  All rights reserved.
*
* Please refer to the NVIDIA end user license agreement (EULA) associated
* with this source code for terms and conditions that govern your use of
* this software. Any use, reproduction, disclosure, or distribution of
* this software and related documentation outside the terms of the EULA
* is strictly prohibited.
*
*/

// Simple 3D volume renderer

#ifndef _VOLUMERENDER_KERNEL_CU_
#define _VOLUMERENDER_KERNEL_CU_

#include <helper_cuda.h>
#include <helper_math.h>

#include "macros.h"
#include "parameters.h"
#include "planckmapping.h"
#include "structures.h"
#include "device_parameters.cuh"
#include "device_functions.cuh"

cudaArray *d_volumeArray = 0;
cudaArray *d_transferFuncArray;
cudaArray *d_visualHullArray = 0;

__constant__ float3x4 c_invViewMatrix;  // inverse view matrix
__constant__ float3 c_nearFar[8];

//texture<VolumeType, 3, cudaReadModeNormalizedFloat> tex;         // 3D texture
texture<VolumeType, 3, cudaReadModeElementType> tex;         // 3D texture
texture<float4, 1, cudaReadModeElementType>         transferTex; // 1D transfer function texture


__device__ static
uint calcIndex(int x, int y, int z, cudaExtent gridSize)
{
	//x = min(x, gridSize.width - 1);
	//y = min(y, gridSize.height - 1);
	//z = min(z, gridSize.depth - 1);
	x = min(x, int(gridSize.width - 1));
	y = min(y, int(gridSize.height - 1));
	z = min(z, int(gridSize.depth - 1));
	x = max(0, x);
	y = max(0, y);
	z = max(0, z);
	uint i = (z*gridSize.width*gridSize.height) + (y*gridSize.width) + x;
	return  i;
}
__device__ static
uint3 calcPos(uint index, uint3 voxelSize)
{
	uint3 voxelPos;
	voxelPos.x = index % voxelSize.x;
	int temp = index / voxelSize.x;
	voxelPos.y = (temp % voxelSize.y);
	voxelPos.z = (temp / voxelSize.y);
	return voxelPos;
}

__device__ 
void d_isotropic_scattering(float& color, float& voxelLength, float3& pos, 
float3& boxMin, float3& boxMax, cudaExtent& extent)
{
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x - voxelLength) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);

	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x + voxelLength) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y - voxelLength) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y + voxelLength) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z - voxelLength) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z + voxelLength) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);

	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x + voxelLength) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y + voxelLength) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z + voxelLength) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x - voxelLength) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y + voxelLength) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z + voxelLength) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x + voxelLength) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y - voxelLength) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z + voxelLength) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x + voxelLength) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y + voxelLength) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z - voxelLength) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x - voxelLength) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y - voxelLength) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z + voxelLength) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x + voxelLength) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y - voxelLength) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z - voxelLength) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x - voxelLength) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y + voxelLength) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z - voxelLength) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
	color += 1.f/ 4 / PI * tex3D(tex,
		(0.5f + ((pos.x - voxelLength) - boxMin.x) / (boxMax.x - boxMin.x) * (extent.width - 1)) / extent.width,
		(0.5f + ((pos.y - voxelLength) - boxMin.y) / (boxMax.y - boxMin.y) * (extent.height - 1)) / extent.height,
		(0.5f + ((pos.z - voxelLength) - boxMin.z) / (boxMax.z - boxMin.z) * (extent.depth - 1)) / extent.depth);
}

__global__ void
d_render_temperature(uint *d_output, uint imageW, uint imageH,
float density, float brightness,
float transferOffset, float transferScale,
int maxSteps, float tstep, float voxelLength,
float3 boxMin, float3 boxMax, cudaExtent volumeSize)
{
	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;

	if ((x >= imageW) || (y >= imageH)) return;

	// calculate eye ray in world space
	Ray eyeRay;
	eyeRay.o = c_nearFar[0] + x * (c_nearFar[3] - c_nearFar[0]) / imageW + y * (c_nearFar[1] - c_nearFar[0]) / imageH;
	eyeRay.d = c_nearFar[4] + x * (c_nearFar[7] - c_nearFar[4]) / imageW + y * (c_nearFar[5] - c_nearFar[4]) / imageH;
	eyeRay.d = normalize(eyeRay.d - eyeRay.o);

	// find intersection with box
	float tnear, tfar;
	int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);

	if (!hit) return;

	// march along ray from front to back, accumulating color
	float4 sum = make_float4(0.0f, 0.f, 0.f, 1.f);
	float t = tnear;
	float3 pos = eyeRay.o + eyeRay.d*tnear;
	float3 step = eyeRay.d*tstep;
	float3 normPos;
	float sample;

	for (int i = 0; i<maxSteps; i++)
	{

		normPos.x = (0.5f + ((pos.x) - boxMin.x) / (boxMax.x - boxMin.x) * (volumeSize.width - 1)) / volumeSize.width;
		normPos.y = (0.5f + ((pos.y) - boxMin.y) / (boxMax.y - boxMin.y) * (volumeSize.height - 1)) / volumeSize.height;
		normPos.z = (0.5f + ((pos.z) - boxMin.z) / (boxMax.z - boxMin.z) * (volumeSize.depth - 1)) / volumeSize.depth;

		// read from 3D texture
		sample = tex3D(tex, normPos.x, normPos.y, normPos.z);
		if (sample > CLAMP_TEMPERATURE_LOW)
		{
			d_isotropic_scattering(sample, voxelLength, pos, boxMin, boxMax, volumeSize);

			float4 col = tex1D(transferTex, 
				(sample - CLAMP_TEMPERATURE_LOW) / (CLAMP_TEMPERATURE_HIGH - CLAMP_TEMPERATURE_LOW));


			// under operator for front-to-back blending
			sum.x = sum.w * col.w * col.x + sum.x;
			sum.y = sum.w * col.w * col.y + sum.y;
			sum.z = sum.w * col.w * col.z + sum.z;
			sum.w = (1 - col.w) * sum.w;

		}

		t += tstep;

		if (t > tfar) break;

		pos += step;
	}


	// write output color
	//sum *= brightness;
	d_output[y*imageW + x] = rgbaFloatToInt(sum);
}



__global__ void
d_render_color(uint *d_output, uint imageW, uint imageH,
float density, float brightness,
float transferOffset, float transferScale,
int maxSteps, float tstep,
float3 boxMin, float3 boxMax, cudaExtent volumeSize, 
float voxelLength, int channel)
{
	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;

	if ((x >= imageW) || (y >= imageH)) return;


	// calculate eye ray in world space
	Ray eyeRay;
	eyeRay.o = c_nearFar[0] + x * (c_nearFar[3] - c_nearFar[0]) / imageW + y * (c_nearFar[1] - c_nearFar[0]) / imageH;
	eyeRay.d = c_nearFar[4] + x * (c_nearFar[7] - c_nearFar[4]) / imageW + y * (c_nearFar[5] - c_nearFar[4]) / imageH;
	eyeRay.d = normalize(eyeRay.d - eyeRay.o);

	// find intersection with box
	float tnear, tfar;
	int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);

	if (!hit) return;

	//if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

	// march along ray from front to back, accumulating color
	float4 sum = make_float4(0.0f, 0.f, 0.f, 1.f);
	float t = tnear;
	float3 pos = eyeRay.o + eyeRay.d*tnear;
	float3 step = eyeRay.d*tstep;
	float3 normPos;
	float sample = 0;

	for (int i = 0; i<maxSteps; i++)
	{
		// read from 3D texture
		// remap position to [0, 1] coordinates
		// +0.5是由于cuda 3dtex的存储方式
		//float sample = tex3D(tex, (0.5 + (pos.x - boxMin.x) / (boxMax.x - boxMin.x) * (volumeSize.width - 1)) / volumeSize.width, 
		//	(0.5 + (pos.y - boxMin.y) / (boxMax.y - boxMin.y) * (volumeSize.height - 1)) / volumeSize.height,
		//	(0.5 + (pos.z - boxMin.z) / (boxMax.z - boxMin.z) * (volumeSize.depth - 1)) / volumeSize.depth);


		normPos.x = (0.5f + ((pos.x) - boxMin.x) / (boxMax.x - boxMin.x) * (volumeSize.width - 1)) / volumeSize.width;
		normPos.y = (0.5f + ((pos.y) - boxMin.y) / (boxMax.y - boxMin.y) * (volumeSize.height - 1)) / volumeSize.height;
		normPos.z = (0.5f + ((pos.z) - boxMin.z) / (boxMax.z - boxMin.z) * (volumeSize.depth - 1)) / volumeSize.depth;

		sample = tex3D(tex, normPos.x, normPos.y, normPos.z);

		if (sample > 0)
		{

			d_isotropic_scattering(sample, voxelLength, pos, boxMin, boxMax, volumeSize);

			// lookup in transfer function texture
			float4 col = { 0.0, 0.0, 0.0, ALPHA };

			if (0 == channel)
			{
				col.x = sample;
			}
			else if (1 == channel)
			{
				col.y = sample;
			}
			else if (2 == channel)
			{
				col.z = sample;
			}

			//col.w *= density;

			// "under" operator for back-to-front blending
			//sum = lerp(sum, col, col.w);
			// sl: maybe "under" operator for front-to-back blending
			//sum = (1 - sum.w) * col + sum;


			//// pre-multiply alpha
			//col.x *= col.w;
			//col.y *= col.w;
			//col.z *= col.w;
			//// "over" operator for front-to-back blending
			//sum = sum + col*(1.0f - sum.w);



			// under operator for front-to-back blending
			sum.x = sum.w * col.w * col.x + sum.x;
			sum.y = sum.w * col.w * col.y + sum.y;
			sum.z = sum.w * col.w * col.z + sum.z;
			sum.w = (1 - col.w) * sum.w;

			//// exit early if opaque
			//if (sum.w > opacityThreshold)
			//	break;

		}

		t += tstep;

		if (t > tfar) break;

		pos += step;
	}

	sum *= brightness;

	// write output color
	d_output[y*imageW + x] = rgbaFloatToInt(sum);
}

__global__ void
d_render_color_temperature(uint *d_output, uint imageW, uint imageH,
		float gMinGreen, float gMaxGreen,
float density, float brightness,
float transferOffset, float transferScale,
int maxSteps, float tstep,
float3 boxMin, float3 boxMax, cudaExtent volumeSize,
float voxelLength, int channel)
{
	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;

	if ((x >= imageW) || (y >= imageH)) return;


	// calculate eye ray in world space
	Ray eyeRay;
	eyeRay.o = c_nearFar[0] + x * (c_nearFar[3] - c_nearFar[0]) / imageW + y * (c_nearFar[1] - c_nearFar[0]) / imageH;
	eyeRay.d = c_nearFar[4] + x * (c_nearFar[7] - c_nearFar[4]) / imageW + y * (c_nearFar[5] - c_nearFar[4]) / imageH;
	eyeRay.d = normalize(eyeRay.d - eyeRay.o);

	// find intersection with box
	float tnear, tfar;
	int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);

	if (!hit) return;

	//if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

	// march along ray from front to back, accumulating color
	float4 sum = make_float4(0.0f, 0.f, 0.f, 1.f);
	float t = tnear;
	float3 pos = eyeRay.o + eyeRay.d*tnear;
	float3 step = eyeRay.d*tstep;
	float3 normPos;
	float sample = 0;

	for (int i = 0; i<maxSteps; i++)
	{
		// read from 3D texture
		// remap position to [0, 1] coordinates
		// +0.5是由于cuda 3dtex的存储方式
		//float sample = tex3D(tex, (0.5 + (pos.x - boxMin.x) / (boxMax.x - boxMin.x) * (volumeSize.width - 1)) / volumeSize.width, 
		//	(0.5 + (pos.y - boxMin.y) / (boxMax.y - boxMin.y) * (volumeSize.height - 1)) / volumeSize.height,
		//	(0.5 + (pos.z - boxMin.z) / (boxMax.z - boxMin.z) * (volumeSize.depth - 1)) / volumeSize.depth);


		normPos.x = (0.5f + ((pos.x) - boxMin.x) / (boxMax.x - boxMin.x) * (volumeSize.width - 1)) / volumeSize.width;
		normPos.y = (0.5f + ((pos.y) - boxMin.y) / (boxMax.y - boxMin.y) * (volumeSize.height - 1)) / volumeSize.height;
		normPos.z = (0.5f + ((pos.z) - boxMin.z) / (boxMax.z - boxMin.z) * (volumeSize.depth - 1)) / volumeSize.depth;

		sample = tex3D(tex, normPos.x, normPos.y, normPos.z);

		if (sample > gMinGreen)
		{

			d_isotropic_scattering(sample, voxelLength, pos, boxMin, boxMax, volumeSize);

			float4 transferCol = tex1D(transferTex, 
				(sample - gMinGreen) / (gMaxGreen - gMinGreen));

			float4 col = { transferCol.x, sample, transferCol.z, transferCol.w};


			// under operator for front-to-back blending
			sum.x = sum.w * col.w * col.x + sum.x;
			sum.y = sum.w * col.w * col.y + sum.y;
			sum.z = sum.w * col.w * col.z + sum.z;
			sum.w = (1 - col.w) * sum.w;
		}

		t += tstep;

		if (t > tfar) break;

		pos += step;
	}

	sum *= brightness;

	// write output color
	d_output[y*imageW + x] = rgbaFloatToInt(sum);
}



extern "C"
void setTextureFilterMode(bool bLinearFilter)
{
	tex.filterMode = bLinearFilter ? cudaFilterModeLinear : cudaFilterModePoint;
}

extern "C"
void initColorTemperatureTex()
{
	//ofstream outStream;
	//outStream.open("data/color_temperature.txt");
	PlanckMapping pm;
	size_t transferFuncSize = CLAMP_TEMPERATURE_HIGH - CLAMP_TEMPERATURE_LOW + 1;
	float4* transferFunc = (float4*)malloc(transferFuncSize * sizeof(float4));
	gMinGreenIntensity = PlanckMapping::bbSpectrum[0].c[1];
	gMaxGreenIntensity = PlanckMapping::bbSpectrum[transferFuncSize - 1].c[1];
	transferFunc[0].x = PlanckMapping::bbSpectrum[0].c[0];
	transferFunc[0].y = CLAMP_TEMPERATURE_LOW;
	transferFunc[0].z = PlanckMapping::bbSpectrum[0].c[2];
	transferFunc[0].w = 0.01f;
	transferFunc[transferFuncSize - 1].x = PlanckMapping::bbSpectrum[transferFuncSize - 1].c[0];
	transferFunc[transferFuncSize - 1].y = CLAMP_TEMPERATURE_HIGH;
	transferFunc[transferFuncSize - 1].z = PlanckMapping::bbSpectrum[transferFuncSize - 1].c[2];
	transferFunc[transferFuncSize - 1].w = 0.05f;
	float step = 1.f * (gMaxGreenIntensity - gMinGreenIntensity) / (CLAMP_TEMPERATURE_HIGH - CLAMP_TEMPERATURE_LOW);
	int currentIndex = 0;
	for (int i = 1; i < transferFuncSize - 1; ++i)
	{
		for (int j = currentIndex; j < transferFuncSize - 1; ++j)
		{
			if (gMinGreenIntensity + i * step <= PlanckMapping::bbSpectrum[j + 1].c[1])
			{
				transferFunc[i].x = (gMinGreenIntensity + i * step - PlanckMapping::bbSpectrum[j].c[1])
					/ (PlanckMapping::bbSpectrum[j + 1].c[1] - PlanckMapping::bbSpectrum[j].c[1])
					* (PlanckMapping::bbSpectrum[j + 1].c[0] - PlanckMapping::bbSpectrum[j].c[0])
					+ PlanckMapping::bbSpectrum[j].c[0];
				transferFunc[i].y = (gMinGreenIntensity + i * step - PlanckMapping::bbSpectrum[j].c[1])
					/ (PlanckMapping::bbSpectrum[j + 1].c[1] - PlanckMapping::bbSpectrum[j].c[1])
					+ CLAMP_TEMPERATURE_LOW + j;
				transferFunc[i].z = (gMinGreenIntensity + i * step - PlanckMapping::bbSpectrum[j].c[1])
					/ (PlanckMapping::bbSpectrum[j + 1].c[1] - PlanckMapping::bbSpectrum[j].c[1])
					* (PlanckMapping::bbSpectrum[j + 1].c[2] - PlanckMapping::bbSpectrum[j].c[2])
					+ PlanckMapping::bbSpectrum[j].c[2];
				transferFunc[i].w = 0.01f;
				break;
			}
			else
			{
				++currentIndex;
			}
		}
		/*outStream << "green = " << gMinGreenIntensity + i * step << "\t"
			<< "red = " << transferFunc[i].x << "\t"
			<< "blue = " << transferFunc[i].z << "\t"
			<< "temperature = " << transferFunc[i].y << endl;*/
	}

	//outStream.close();
	//getchar();
	cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<float4>();
	cudaArray *d_transferFuncArray;
	checkCudaErrors(cudaMallocArray(&d_transferFuncArray, &channelDesc2, transferFuncSize, 1));
	checkCudaErrors(cudaMemcpyToArray(d_transferFuncArray, 0, 0, transferFunc, transferFuncSize * sizeof(float4), cudaMemcpyHostToDevice));
	//checkCudaErrors(cudaMallocArray(&d_transferFuncArray, &channelDesc2, sizeof(transferFunc) / sizeof(float4), 1));
	//checkCudaErrors(cudaMemcpyToArray(d_transferFuncArray, 0, 0, transferFunc, sizeof(transferFunc), cudaMemcpyHostToDevice));

	transferTex.filterMode = cudaFilterModeLinear;
	transferTex.normalized = true;    // access with normalized texture coordinates
	transferTex.addressMode[0] = cudaAddressModeClamp;   // wrap texture coordinates

	// Bind the array to the texture
	checkCudaErrors(cudaBindTextureToArray(transferTex, d_transferFuncArray, channelDesc2));
}


extern "C"
void initTransferTex()
{
	size_t transferFuncSize = CLAMP_TEMPERATURE_HIGH - CLAMP_TEMPERATURE_LOW + 1;
	float4* transferFunc = (float4*)malloc(transferFuncSize * sizeof(float4));
	for (int i = 0; i < transferFuncSize; ++i)
	{
		transferFunc[i].x = PlanckMapping::bbSpectrum[i].c[0] ;
		transferFunc[i].y = PlanckMapping::bbSpectrum[i].c[1];
		transferFunc[i].z = PlanckMapping::bbSpectrum[i].c[2];
		// to do
		transferFunc[i].w = 0.05f;

	}

	cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<float4>();
	cudaArray *d_transferFuncArray;
	checkCudaErrors(cudaMallocArray(&d_transferFuncArray, &channelDesc2, transferFuncSize, 1));
	checkCudaErrors(cudaMemcpyToArray(d_transferFuncArray, 0, 0, transferFunc, transferFuncSize * sizeof(float4), cudaMemcpyHostToDevice));
	//checkCudaErrors(cudaMallocArray(&d_transferFuncArray, &channelDesc2, sizeof(transferFunc) / sizeof(float4), 1));
	//checkCudaErrors(cudaMemcpyToArray(d_transferFuncArray, 0, 0, transferFunc, sizeof(transferFunc), cudaMemcpyHostToDevice));

	transferTex.filterMode = cudaFilterModeLinear;
	transferTex.normalized = true;    // access with normalized texture coordinates
	transferTex.addressMode[0] = cudaAddressModeClamp;   // wrap texture coordinates

	// Bind the array to the texture
	checkCudaErrors(cudaBindTextureToArray(transferTex, d_transferFuncArray, channelDesc2));
}

extern "C"
void initCuda(void *h_volume, cudaExtent volumeSize)
{
	// create 3D array
	cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<VolumeType>();
	checkCudaErrors(cudaMalloc3DArray(&d_volumeArray, &channelDesc, volumeSize));

	// copy data to 3D array
	cudaMemcpy3DParms copyParams = { 0 };
	copyParams.srcPtr = make_cudaPitchedPtr(h_volume, volumeSize.width*sizeof(VolumeType), volumeSize.width, volumeSize.height);
	copyParams.dstArray = d_volumeArray;
	copyParams.extent = volumeSize;
	copyParams.kind = cudaMemcpyHostToDevice;
	checkCudaErrors(cudaMemcpy3D(&copyParams));

	// set texture parameters
	tex.normalized = true;                      // access with normalized texture coordinates
	//tex.filterMode = cudaFilterModePoint;      // linear interpolation only support float type
	tex.filterMode = cudaFilterModeLinear;      // linear interpolation only support float type
	tex.addressMode[0] = cudaAddressModeClamp;  // clamp texture coordinates
	tex.addressMode[1] = cudaAddressModeClamp;

	// bind array to 3D texture
	checkCudaErrors(cudaBindTextureToArray(tex, d_volumeArray, channelDesc));

	// create transfer function texture
	//float4 transferFunc[] =
	//{
	//	{ 0.0, 0.0, 0.0, 0.0, },
	//	{ 1.0, 0.0, 0.0, 1.0, },
	//	{ 1.0, 0.5, 0.0, 1.0, },
	//	{ 1.0, 1.0, 0.0, 1.0, },
	//	{ 0.0, 1.0, 0.0, 1.0, },
	//	{ 0.0, 1.0, 1.0, 1.0, },
	//	{ 0.0, 0.0, 1.0, 1.0, },
	//	{ 1.0, 0.0, 1.0, 1.0, },
	//	{ 0.0, 0.0, 0.0, 0.0, },
	//};
	size_t transferFuncSize = CLAMP_TEMPERATURE_HIGH - CLAMP_TEMPERATURE_LOW + 1;
	float4* transferFunc = (float4*)malloc(transferFuncSize * sizeof(float4));
	for (int i = 0; i < transferFuncSize; ++i)
	{
		transferFunc[i].x = PlanckMapping::bbSpectrum[i].c[0];
		transferFunc[i].y = PlanckMapping::bbSpectrum[i].c[1];
		transferFunc[i].z = PlanckMapping::bbSpectrum[i].c[2];
		// to do
		transferFunc[i].w = 0.2f;

	}

	cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<float4>();
	cudaArray *d_transferFuncArray;
	checkCudaErrors(cudaMallocArray(&d_transferFuncArray, &channelDesc2,transferFuncSize, 1));
	checkCudaErrors(cudaMemcpyToArray(d_transferFuncArray, 0, 0, transferFunc, transferFuncSize * sizeof(float4), cudaMemcpyHostToDevice));
	//checkCudaErrors(cudaMallocArray(&d_transferFuncArray, &channelDesc2, sizeof(transferFunc) / sizeof(float4), 1));
	//checkCudaErrors(cudaMemcpyToArray(d_transferFuncArray, 0, 0, transferFunc, sizeof(transferFunc), cudaMemcpyHostToDevice));

	transferTex.filterMode = cudaFilterModeLinear;
	transferTex.normalized = true;    // access with normalized texture coordinates
	transferTex.addressMode[0] = cudaAddressModeClamp;   // wrap texture coordinates

	// Bind the array to the texture
	checkCudaErrors(cudaBindTextureToArray(transferTex, d_transferFuncArray, channelDesc2));
}

extern "C"
void bindVisualHullTexture(void *h_volume, cudaExtent volumeSize)
{
	// create 3D array
	cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<unsigned short>();
	checkCudaErrors(cudaMalloc3DArray(&d_visualHullArray, &channelDesc, volumeSize));

	// copy data to 3D array
	cudaMemcpy3DParms copyParams = { 0 };
	copyParams.srcPtr = make_cudaPitchedPtr(h_volume, volumeSize.width*sizeof(unsigned short), volumeSize.width, volumeSize.height);
	copyParams.dstArray = d_visualHullArray;
	copyParams.extent = volumeSize;
	copyParams.kind = cudaMemcpyHostToDevice;
	checkCudaErrors(cudaMemcpy3D(&copyParams));

	// set texture parameters
	tex.normalized = true;                      // access with normalized texture coordinates
	tex.filterMode = cudaFilterModeLinear;      // linear interpolation
	tex.addressMode[0] = cudaAddressModeClamp;  // clamp texture coordinates
	tex.addressMode[1] = cudaAddressModeClamp;

	// bind array to 3D texture
	checkCudaErrors(cudaBindTextureToArray(tex, d_visualHullArray, channelDesc));

}
extern "C"
void bindVolumeTexture(void *h_volume, cudaExtent volumeSize)
{
	// create 3D array
	cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<VolumeType>();
	checkCudaErrors(cudaMalloc3DArray(&d_volumeArray, &channelDesc, volumeSize));

	// copy data to 3D array
	cudaMemcpy3DParms copyParams = { 0 };
	copyParams.srcPtr = make_cudaPitchedPtr(h_volume, volumeSize.width*sizeof(VolumeType), volumeSize.width, volumeSize.height);
	copyParams.dstArray = d_volumeArray;
	copyParams.extent = volumeSize;
	copyParams.kind = cudaMemcpyHostToDevice;
	checkCudaErrors(cudaMemcpy3D(&copyParams));

	// set texture parameters
	tex.normalized = true;                      // access with normalized texture coordinates
	tex.filterMode = cudaFilterModeLinear;      // linear interpolation
	tex.addressMode[0] = cudaAddressModeClamp;  // clamp texture coordinates
	tex.addressMode[1] = cudaAddressModeClamp;

	// bind array to 3D texture
	checkCudaErrors(cudaBindTextureToArray(tex, d_volumeArray, channelDesc));

}

extern "C"
void upgradeVolumeData(cudaExtent volumeExtent)
{
	checkCudaErrors(cudaFreeArray(d_volumeArray));
	checkCudaErrors(cudaMemcpy(gVolumeColor, d_volume_color,
			volumeExtent.width*sizeof(VolumeType)* volumeExtent.width * volumeExtent.height, cudaMemcpyDeviceToHost));
	bindVolumeTexture(gVolumeColor, volumeExtent);
}

extern "C"
void freeCudaBuffers()
{
	if (0 != d_visualHullArray)
	{
		checkCudaErrors(cudaFreeArray(d_visualHullArray));
	}
	checkCudaErrors(cudaFreeArray(d_volumeArray));
	checkCudaErrors(cudaFreeArray(d_transferFuncArray));
}


extern "C"
void render_temperature_kernel(dim3 gridSize, dim3 blockSize, uint *d_output, uint imageW, uint imageH,
float density, float brightness, float transferOffset, float transferScale, 
int samplingMaxSteps, float samplingStep,
float3 boxMin, float3 boxMax, cudaExtent volumeSize)
{
	d_render_temperature << <gridSize, blockSize >> >(d_output, imageW, imageH, density,
		brightness, transferOffset, transferScale, 
		samplingMaxSteps, samplingStep, gVoxelLength,
		boxMin, boxMax, volumeSize);
}


extern "C"
void render_color_kernel(dim3 gridSize, dim3 blockSize, uint *d_output, uint imageW, uint imageH,
float density, float brightness, float transferOffset, float transferScale, 
int samplingMaxSteps, float samplingStep,
float3 boxMin, float3 boxMax, cudaExtent volumeSize, int channel)
{ 
	d_render_color<< <gridSize, blockSize >> >(d_output, imageW, imageH, density,
		brightness, transferOffset, transferScale, 
		samplingMaxSteps, samplingStep,
		boxMin, boxMax, volumeSize, gVoxelLength, gChannel);
	cudaThreadSynchronize();
}


extern "C"
void render_color_temperature_kernel(dim3 gridSize, dim3 blockSize, uint *d_output, uint imageW, uint imageH,
float density, float brightness, float transferOffset, float transferScale,
int samplingMaxSteps, float samplingStep,
float3 boxMin, float3 boxMax, cudaExtent volumeSize, int channel)
{
	d_render_color_temperature << <gridSize, blockSize >> >(d_output, imageW, imageH, 
			gMinGreenIntensity, gMaxGreenIntensity,
			density,brightness, transferOffset, transferScale,
		samplingMaxSteps, samplingStep,
		boxMin, boxMax, volumeSize, gVoxelLength, gChannel);
	cudaThreadSynchronize();
}

extern "C"
void copyInvViewMatrix(float *invViewMatrix, size_t sizeofMatrix)
{
	checkCudaErrors(cudaMemcpyToSymbol(c_invViewMatrix, invViewMatrix, sizeofMatrix));
}

extern "C"
void copyNearFar(float* nearFar, size_t sizeofNearFar)
{
	checkCudaErrors(cudaMemcpyToSymbol(c_nearFar, nearFar, sizeofNearFar));
}

__device__ void d_adjust_temperature(float* volumeData, uchar* visualHull, float3 normPos, cudaExtent voxelSize,
	int xoffset, int yoffset, int zoffset, float midDis,
	float residual, float transparency, int sampleInd, float laplace, float src,
	int numberCamera, unsigned long iteration)
{
	float3 dstPos, dstNormPos;

	uint idx;
	// 8 neighbor key points
	// key point relative coordinate
	dstPos.x = int(normPos.x * (voxelSize.width) - 0.5f) + xoffset;
	dstPos.y = int(normPos.y * (voxelSize.height) - 0.5f) + yoffset;
	dstPos.z = int(normPos.z * (voxelSize.depth) - 0.5f) + zoffset;


	idx = calcIndex(
		int(dstPos.x),
		int(dstPos.y),
		int(dstPos.z),
		voxelSize);


	//if (int(visualHull[idx]) == numberCamera)
	{

		// key point 3d tex normalized coordinate
		dstNormPos.x = (dstPos.x + 0.5f) / voxelSize.width;
		dstNormPos.y = (dstPos.y + 0.5f) / voxelSize.height;
		dstNormPos.z = (dstPos.z + 0.5f) / voxelSize.depth;



		float currentColor = volumeData[idx];

		float dis = dot(dstNormPos - normPos,
			dstNormPos - normPos);

		//printf("midDis = %f\n", midDis);
		//float delta = 1000* 0.001f * residual * transparency / exp(dis);
		//float delta = 20 * 0.1f * 0.01f * residual * pow(1 - transparency, sampleInd - 1) / exp(dis);
		//float delta = 10 * 20 * 0.1f * (0.01f * residual + 0.00f * laplace) * transparency / exp(dis);
		//float delta = 200 * 0.1f * 0.01f * residual / (1.f + exp(40.f * midDis - 3.f)) / exp(dis);
		//float delta = 5 * 0.1 * 2 / (1.f + exp(iteration - 6.f)) * residual / (1.f + exp(40.f * midDis - 3.f)) / exp(dis);
		//float delta = 1 * 0.1f * 0.01f * residual * residual / (1.f + exp(40.f * midDis - 3.f)) / exp(dis);
		//printf("logistic = %f\n", 1.f / (1.f + exp(10.f * midDis - 5.f)));
		//float delta =  2* 0.1f * 0.01f * residual / exp(dis);
		float delta = 20 * 0.1f * 0.01f * residual;

		//float delta;
		//curandState s;
		//curand_init(idx, int(iterTime * transparency), 0, &s);
		//SDK部分delta值
		//float delta = 100 * 0.5 */*(curand_uniform(&s) + 1) * */(0.001*laplace + 0.001f * residual)
		//		* exp((1 - (transparency - 0.6)*(transparency - 0.6)) + 0.64 - 100.0f * dis /*- 0.1f * iterTime*/);
		//float delta = residual * 10.1;
		//if (fabs(residual) > 10e-3)
		//{
		//	delta = 10000 * 0.05 */*(curand_uniform(&s) + 1) * */(0.001*laplace + 0.001f * residual)
		//		* exp((1 - (transparency - 0.6)*(transparency - 0.6)) + 0.64 - 100.0f * dis /*- 0.1f * iterTime*/);
		//}
		//else
		//{
		//	delta = -1.f;
		//}


		//volumeData[idx] += delta;
		currentColor = currentColor + delta;
		//volumeData[idx] = currentColor + delta;
		//atomicAdd(volumeData + idx, delta);
		//// may be improved
		if (currentColor >= 300 && currentColor <= CLAMP_TEMPERATURE_HIGH)
		{
			//atomicAdd(volumeData + idx, delta);
			volumeData[idx] += delta;
		}
		else 
		if (currentColor < 300)
		{
			volumeData[idx] = 300;
		}
		else 
		{
			volumeData[idx] = CLAMP_TEMPERATURE_HIGH;
		}
	}
}


__device__ void d_adjust_color(float* volumeData, uchar* visualHull, float3 normPos, cudaExtent voxelSize,
	int xoffset, int yoffset, int zoffset,
	float residual, float transparency, int sampleInd, float laplace, float src,
	int numberCamera)
{
	float3 dstPos, dstNormPos;

	uint idx;
	// 8 neighbor key points
	// key point relative coordinate
	dstPos.x = int(normPos.x * (voxelSize.width) - 0.5f) + xoffset;
	dstPos.y = int(normPos.y * (voxelSize.height) - 0.5f) + yoffset;
	dstPos.z = int(normPos.z * (voxelSize.depth) - 0.5f) + zoffset;


	idx = calcIndex(
		int(dstPos.x),
		int(dstPos.y),
		int(dstPos.z),
		voxelSize);
	//int tag = int(visualHull[idx]);
#if ENABLE_ONE_VIEW_TEST_COLOR
	if (int(visualHull[idx]) == numberCamera - 1)
#else
	if (int(visualHull[idx]) == numberCamera)
#endif
	//if (tag == numberCamera)
	{

		// key point 3d tex normalized coordinate
		dstNormPos.x = (dstPos.x + 0.5f) / voxelSize.width;
		dstNormPos.y = (dstPos.y + 0.5f) / voxelSize.height;
		dstNormPos.z = (dstPos.z + 0.5f) / voxelSize.depth;



		float currentColor = volumeData[idx];

		float dis = dot(dstNormPos - normPos,
			dstNormPos - normPos);


		//float delta = (0.001f * residual + 0.00f * laplace) * transparency / exp(dis);
		float delta = 0.001f * residual * transparency / exp(dis);


		//volumeData[idx] += delta;
		currentColor = currentColor + delta;
		//volumeData[idx] = currentColor + delta;
		//atomicAdd(volumeData + idx, delta);
		//// may be improved
		if (currentColor >= 0 && currentColor <= 255)
		{
			//atomicAdd(volumeData + idx, delta);
			volumeData[idx] += delta;
		}
		else if (currentColor < 0)
		{
			volumeData[idx] = 0;
		}
		else
		{
			volumeData[idx] = 255;
		}
	}
}



__global__ void
d_reconstruct_color(float* volumeData, uchar* visualHull, 
int numberCamera, int cameraIndex, int widthStep, int nChannels,
char* imageData, uint *d_output,  uint imageW, uint imageH,
float density, float brightness,
float transferOffset, float transferScale,
int maxSteps, float tstep,
float3 boxMin, float3 boxMax,
cudaExtent volumeSize, int channel)
{

	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;

	if ((x >= imageW) || (y >= imageH)) return;

	float4 color = IntToRgbaFloat(d_output[y * imageW + x]) ;
	float cur = 0.f;
	if (0 == channel)
	{
		cur = color.x;
	}
	else if (1 == channel)
	{
		cur = color.y;
	}
	else if (2 == channel)
	{
		cur = color.z;
	}

	float src = uchar(imageData[cameraIndex * imageH * widthStep 
		+ (imageH - y - 1) * widthStep + x * nChannels + 2 - channel]);

	float residual = src - cur;
	float lap = 0.f;
	if (x < imageW - 1 && x > 0 && y < imageH - 1 && y > 0)
	{
		float curUp = 0.f;
		float curDown = 0.f;
		float curLeft = 0.f;
		float curRight = 0.f;
		if (0 == channel)
		{
			curUp		= IntToRgbaFloat(d_output[(y - 1) * imageW + x]).x;
			curDown		= IntToRgbaFloat(d_output[(y + 1) * imageW + x]).x;
			curLeft		= IntToRgbaFloat(d_output[y * imageW + x - 1]).x;
			curRight	= IntToRgbaFloat(d_output[y * imageW + x + 1]).x;
		}
		else if (1 == channel)
		{
			curUp		= IntToRgbaFloat(d_output[(y - 1) * imageW + x]).y;
			curDown		= IntToRgbaFloat(d_output[(y + 1) * imageW + x]).y;
			curLeft		= IntToRgbaFloat(d_output[y * imageW + x - 1]).y;
			curRight	= IntToRgbaFloat(d_output[y * imageW + x + 1]).y;
		}
		else if (2 == channel)
		{
			curUp		= IntToRgbaFloat(d_output[(y - 1) * imageW + x]).z;
			curDown		= IntToRgbaFloat(d_output[(y + 1) * imageW + x]).z;
			curLeft		= IntToRgbaFloat(d_output[y * imageW + x - 1]).z;
			curRight	= IntToRgbaFloat(d_output[y * imageW + x + 1]).z;
		}

		float srcUp = uchar(imageData[cameraIndex * imageH * widthStep
			+ (imageH - y) * widthStep + x * nChannels + 2 - channel]);
		float srcDown = uchar(imageData[cameraIndex * imageH * widthStep
			+ (imageH - y - 2) * widthStep + x * nChannels + 2 - channel]);
		float srcLeft = uchar(imageData[cameraIndex * imageH * widthStep
			+ (imageH - y - 1) * widthStep + (x - 1) * nChannels + 2 - channel]);
		float srcRight = uchar(imageData[cameraIndex * imageH * widthStep
			+ (imageH - y - 1) * widthStep + (x + 1) * nChannels + 2 - channel]);

		lap = (srcUp + srcDown + srcLeft + srcRight - 4 * src) - (curUp + curDown + curLeft + curRight - 4 * cur);
	}

	if (fabsf(residual) < 10e-3)
	{
		return;
	}
	//if ((residual < 0 && residual > -10e-3) || (residual > 0 && residual < 10e-3))
	//{
	//	return;
	//}



	// 计算delat值完毕，后续进行强度调整
	// calculate eye ray in world space
	Ray eyeRay;
	eyeRay.o = c_nearFar[0] + x * (c_nearFar[3] - c_nearFar[0]) / imageW + y * (c_nearFar[1] - c_nearFar[0]) / imageH;
	eyeRay.d = c_nearFar[4] + x * (c_nearFar[7] - c_nearFar[4]) / imageW + y * (c_nearFar[5] - c_nearFar[4]) / imageH;
	eyeRay.d = normalize(eyeRay.d - eyeRay.o);

	// find intersection with box
	float tnear, tfar;
	int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);

	if (!hit) return;

	//if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

	// march along ray from front to back, accumulating color
	float t = tnear;
	float3 pos = eyeRay.o + eyeRay.d*tnear;
	float3 step = eyeRay.d*tstep;
	float3 normPos;
	float alpha = 1.f;
	


	for (int i = 0; i<maxSteps; i++)
	{


		// read from 3D texture
		// remap position to [0, 1] coordinates
		// +0.5是由于cuda 3dtex的存储方式

		// 3d tex normalized coordinate
		normPos.x = (0.5f + ((pos.x) - boxMin.x) / (boxMax.x - boxMin.x) * (volumeSize.width - 1)) / volumeSize.width;
		normPos.y = (0.5f + ((pos.y) - boxMin.y) / (boxMax.y - boxMin.y) * (volumeSize.height - 1)) / volumeSize.height;
		normPos.z = (0.5f + ((pos.z) - boxMin.z) / (boxMax.z - boxMin.z) * (volumeSize.depth - 1)) / volumeSize.depth;





		// splatting
		//int scale = myRand() % 3 + 1;
		int scale = 1;
		for (int xoffset = 0; xoffset < scale; ++xoffset)
		{
			for (int yoffset = 0; yoffset < scale; ++yoffset)
			{
				for (int zoffset = 0; zoffset < scale; ++zoffset)
				{
					d_adjust_color(volumeData, visualHull, normPos, volumeSize,
						xoffset, yoffset, zoffset,
						residual, alpha, i, lap, src, numberCamera);
				}
			}
		}

		// "over" operator for front-to-back blending
		//alpha = alpha + density * ALPHA * (1.f - alpha);

		// "under" operator for front-to-back blending
		alpha = (1 - ALPHA) * alpha;


		t += tstep;

		if (t > tfar) break;

		pos += step;
	}

}
__global__ void
d_reconstruct_color_temperature(float* volumeData, uchar* visualHull,
int numberCamera, int cameraIndex, int widthStep, int nChannels,
char* imageData, uint *d_output, uint imageW, uint imageH,
float density, float brightness,
float transferOffset, float transferScale,
int maxSteps, float tstep,
float3 boxMin, float3 boxMax,
cudaExtent volumeSize, int channel)
{

	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;

	if ((x >= imageW) || (y >= imageH)) return;

	float4 color = IntToRgbaFloat(d_output[y * imageW + x]);
	float cur = color.y;

	float src = uchar(imageData[cameraIndex * imageH * widthStep
		+ (imageH - y - 1) * widthStep + x * nChannels + 1]);

	float residual = src - cur;
	float lap = 0.f;
	if (fabsf(residual) < 10e-3)
	{
		return;
	}
	//printf("src = %f, cur = %f\n", src, cur);

	// 计算delat值完毕，后续进行强度调整
	// calculate eye ray in world space
	Ray eyeRay;
	eyeRay.o = c_nearFar[0] + x * (c_nearFar[3] - c_nearFar[0]) / imageW + y * (c_nearFar[1] - c_nearFar[0]) / imageH;
	eyeRay.d = c_nearFar[4] + x * (c_nearFar[7] - c_nearFar[4]) / imageW + y * (c_nearFar[5] - c_nearFar[4]) / imageH;
	eyeRay.d = normalize(eyeRay.d - eyeRay.o);

	// find intersection with box
	float tnear, tfar;
	int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);

	if (!hit) return;


	// march along ray from front to back, accumulating color
	float t = tnear;
	float3 pos = eyeRay.o + eyeRay.d*tnear;
	float3 step = eyeRay.d*tstep;
	float3 normPos;
	float alpha = 1.f;


	for (int i = 0; i<maxSteps; i++)
	{

		// read from 3D texture
		// remap position to [0, 1] coordinates
		// +0.5是由于cuda 3dtex的存储方式

		// 3d tex normalized coordinate
		normPos.x = (0.5f + ((pos.x) - boxMin.x) / (boxMax.x - boxMin.x) * (volumeSize.width - 1)) / volumeSize.width;
		normPos.y = (0.5f + ((pos.y) - boxMin.y) / (boxMax.y - boxMin.y) * (volumeSize.height - 1)) / volumeSize.height;
		normPos.z = (0.5f + ((pos.z) - boxMin.z) / (boxMax.z - boxMin.z) * (volumeSize.depth - 1)) / volumeSize.depth;


		// splatting
		//int scale = myRand() % 3 + 1;
		int scale = 1;
		for (int xoffset = 0; xoffset < scale; ++xoffset)
		{
			for (int yoffset = 0; yoffset < scale; ++yoffset)
			{
				for (int zoffset = 0; zoffset < scale; ++zoffset)
				{
					d_adjust_color(volumeData, visualHull, normPos, volumeSize,
						xoffset, yoffset, zoffset,
						residual, alpha, i, lap, src, numberCamera);
				}
			}
		}

		// "over" operator for front-to-back blending
		//alpha = alpha + density * ALPHA * (1.f - alpha);

		// "under" operator for front-to-back blending
		alpha = (1 - ALPHA) * alpha;

		t += tstep;

		if (t > tfar) break;

		pos += step;
	}

}

__global__ void
d_generate_visual_hull_color_temperature_resize_octree(VisualHullType* visualHull,
	int numberCamera, int widthStep, int nChannels,
	ResizeImage* gResizeImageData, uint imageW, uint imageH,
	cudaExtent volumeSize, view_box visualHullBox, extrinsic_parameters* extrinsic_par, Octree* oct, uint oct_index)
{
	uint len = blockIdx.x * blockDim.x + threadIdx.x;
	if (len >= oct_index) return;

	double _x = visualHullBox.x / volumeSize.width;
	double _y = visualHullBox.y / volumeSize.height;
	double _z = visualHullBox.z / volumeSize.depth;

	double a, b, c;

	a = visualHullBox.o[0] + _x * oct[len].min[0] + _x * oct[len].len / 2;
	b = visualHullBox.o[1] + _y * oct[len].min[1] + _y * oct[len].len / 2;
	c = visualHullBox.o[2] + _z * oct[len].min[2] + _z * oct[len].len / 2;
	//printf("a = %lf, b = %lf, c = %lf\n", a, b, c);


	unsigned short tag = 2;
	double xx, yy;
	int wid, hei;
	int ans;

	for (int i = 0; i < numberCamera; i++) {
		if (!tag) break;

		extrinsic_par[i].d_get_result(&xx, &yy, a, b, c);

		wid = (int)(xx / widthStep*gResizeImageData[i].width);
		hei = (int)(yy / imageH*gResizeImageData[i].height);

		ans = (uchar)gResizeImageData[i].imageData[hei*gResizeImageData[i].width + wid];

		if (volumeSize.width / oct[len].len <= 8) {
			tag = 1;
		}
		else if (ans < 5) {
			tag = 0;
		}
		else if(tag == 2 && ans > 248){
			tag = 2;
		}
		else {
			tag = 1;
		}
	}

	if (tag == 0 || tag == 2 || oct[len].len <= 1) {
		if (tag == 0) {
			oct[len].type = 0;
		}
		else {
			oct[len].type = 2;
		}
	}
	else {
		oct[len].type = 1;
	}

}

__global__ void
d_generate_visual_hull_color_temperature3(VisualHullType* visualHull,
	int numberCamera, int widthStep, int nChannels,
	int* integrogram, uint imageW, uint imageH,
	cudaExtent volumeSize, view_box visualHullBox, extrinsic_parameters* extrinsic_par, Octree* oct, Octree* oct_t,uint* quere, uint oct_index)
{
	uint len = blockIdx.x * blockDim.x + threadIdx.x;
	if (len >= oct_index) return;

	double _x = visualHullBox.x / volumeSize.width;
	double _y = visualHullBox.y / volumeSize.height;
	double _z = visualHullBox.z / volumeSize.depth;

	double a[2], b[2], c[2];

	double f[8][3];
	a[0] = visualHullBox.o[0] + _x * oct[len].min[0];
	a[1] = a[0] + _x * oct[len].len;
	b[0] = visualHullBox.o[1] + _y * oct[len].min[1];
	b[1] = b[0] + _y * oct[len].len;
	c[0] = visualHullBox.o[2] + _z * oct[len].min[2];
	c[1] = c[0] + _z * oct[len].len;
	int zz = 0;
	for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 2; j++) {
			for (int k = 0; k < 2; k++) {
				f[zz][0] = a[i];
				f[zz][1] = b[j];
				f[zz][2] = c[k];
				zz++;
			}
		}
	}

	unsigned short tag = 2;
	double minx, miny, maxx, maxy;
	int iminx, iminy, imaxx, imaxy;
	int ans;

	for (int i = 0; i < numberCamera; i++) {
		if (!tag) break;

		extrinsic_par[i].g_get_M_x_y(&minx, &maxx, &miny, &maxy, f);

		iminx = (int)minx;
		iminy = (int)miny;
		imaxx = (int)maxx + 1;
		imaxy = (int)maxy + 1;

		if (iminx < 0) iminx = 0;
		if (imaxx >= widthStep) imaxx = widthStep - 1;
		if (iminy < 0) iminy = 0;
		if (imaxy >= imageH) imaxy = imageH - 1;

		ans = integrogram[i * imageH * widthStep + iminy * widthStep + iminx] + integrogram[i * imageH * widthStep + imaxy * widthStep + imaxx] - integrogram[i * imageH * widthStep + imaxy * widthStep + iminx] - integrogram[i * imageH * widthStep + iminy * widthStep + imaxx];

		if (ans == 0) {
			tag = 0;
		}
		else if (tag == 2 && ans == (imaxy - iminy)*(imaxx - iminx)) {
			tag = 2;
		}
		else {
			tag = 1;
		}
	}

	if (tag == 0 || tag == 2 || oct[len].len <= 1) {
		if (tag == 0)
			oct[len].type = 0;
		else
			oct[len].type = 2;
	}
	else {
		uint oct_t_index = len * 8;
		quere[oct_t_index] = 1;
		oct_t[oct_t_index++].create(oct[len].min[0], oct[len].min[1], oct[len].min[2], oct[len].len / 2);
		quere[oct_t_index] = 1;
		oct_t[oct_t_index++].create(oct[len].min[0] + oct[len].len / 2, oct[len].min[1], oct[len].min[2], oct[len].len / 2);
		quere[oct_t_index] = 1;
		oct_t[oct_t_index++].create(oct[len].min[0], oct[len].min[1] + oct[len].len / 2, oct[len].min[2], oct[len].len / 2);
		quere[oct_t_index] = 1;
		oct_t[oct_t_index++].create(oct[len].min[0], oct[len].min[1], oct[len].min[2] + oct[len].len / 2, oct[len].len / 2);
		quere[oct_t_index] = 1;
		oct_t[oct_t_index++].create(oct[len].min[0] + oct[len].len / 2, oct[len].min[1] + oct[len].len / 2, oct[len].min[2], oct[len].len / 2);
		quere[oct_t_index] = 1;
		oct_t[oct_t_index++].create(oct[len].min[0] + oct[len].len / 2, oct[len].min[1], oct[len].min[2] + oct[len].len / 2, oct[len].len / 2);
		quere[oct_t_index] = 1;
		oct_t[oct_t_index++].create(oct[len].min[0], oct[len].min[1] + oct[len].len / 2, oct[len].min[2] + oct[len].len / 2, oct[len].len / 2);
		quere[oct_t_index] = 1;
		oct_t[oct_t_index++].create(oct[len].min[0] + oct[len].len / 2, oct[len].min[1] + oct[len].len / 2, oct[len].min[2] + oct[len].len / 2, oct[len].len / 2);
	}

}

// compact voxel array
__global__ void
d_compactOctree(Octree *OctreeArray,Octree *Octrees, uint *voxelOccupied, uint *voxelOccupiedScan, uint numVoxels)
{
	uint blockId = __mul24(blockIdx.y, gridDim.x) + blockIdx.x;
	uint i = __mul24(blockId, blockDim.x) + threadIdx.x;

	if (voxelOccupied[i] && (i < numVoxels))
	{
		OctreeArray[voxelOccupiedScan[i]].create(Octrees[i]);
	}
}

extern "C" void
launch_compactOctree(dim3 grid, dim3 threads, Octree *OctreeArray, Octree *Octrees, uint *voxelOccupied, uint *voxelOccupiedScan, uint numVoxels)
{
	d_compactOctree << <grid, threads >> >(OctreeArray,Octrees, voxelOccupied,
		voxelOccupiedScan, numVoxels);
	//getLastCudaError("compactVoxels failed");
}

__global__ void
d_generate_visual_hull_color_temperature2(VisualHullType* visualHull,
	int numberCamera, int widthStep, int nChannels,
	int* integrogram, uint imageW, uint imageH,
	cudaExtent volumeSize, view_box visualHullBox, extrinsic_parameters* extrinsic_par, Octree* oct, uint oct_index)
{
	uint len = blockIdx.x * blockDim.x + threadIdx.x;
	if (len >= oct_index) return;

	double _x = visualHullBox.x / volumeSize.width;
	double _y = visualHullBox.y / volumeSize.height;
	double _z = visualHullBox.z / volumeSize.depth;

	double a[2], b[2], c[2];

	double f[8][3];
	a[0] = visualHullBox.o[0] + _x * oct[len].min[0];
	a[1] = a[0] + _x * oct[len].len;
	b[0] = visualHullBox.o[1] + _y * oct[len].min[1];
	b[1] = b[0] + _y * oct[len].len;
	c[0] = visualHullBox.o[2] + _z * oct[len].min[2];
	c[1] = c[0] + _z * oct[len].len;
	int zz = 0;
	for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 2; j++) {
			for (int k = 0; k < 2; k++) {
				f[zz][0] = a[i];
				f[zz][1] = b[j];
				f[zz][2] = c[k];
				zz++;
			}
		}
	}

	unsigned short tag = 2;
	double minx, miny, maxx, maxy;
	int iminx, iminy, imaxx, imaxy;
	int ans;

	for (int i = 0; i < numberCamera; i++) {
		if (!tag) break;

		extrinsic_par[i].g_get_M_x_y(&minx, &maxx, &miny, &maxy, f);

		iminx = (int)minx;
		iminy = (int)miny;
		imaxx = (int)maxx + 1;
		imaxy = (int)maxy + 1;

		if (iminx < 0) iminx = 0;
		if (imaxx >= widthStep) imaxx = widthStep - 1;
		if (iminy < 0) iminy = 0;
		if (imaxy >= imageH) imaxy = imageH - 1;

		ans = integrogram[i * imageH * widthStep + iminy * widthStep + iminx] + integrogram[i * imageH * widthStep + imaxy * widthStep + imaxx] - integrogram[i * imageH * widthStep + imaxy * widthStep + iminx] - integrogram[i * imageH * widthStep + iminy * widthStep + imaxx];

		if (ans == 0) {
			tag = 0;
		}
		else if (tag == 2 && ans == (imaxy - iminy)*(imaxx - iminx)) {
			tag = 2;
		}
		else {
			tag = 1;
		}
	}

	if (tag == 0 || tag == 2 || oct[len].len <= 1) {
		if(tag == 0)
			oct[len].type = 0;
		else
			oct[len].type = 2;
	}
	else {
		oct[len].type = 1;
	}

}

__global__ void
d_generate_visual_hull_color_temperature(VisualHullType* visualHull,
int numberCamera, int widthStep, int nChannels,
int* integrogram, uint imageW, uint imageH,
cudaExtent volumeSize, view_box visualHullBox, extrinsic_parameters* extrinsic_par)
{
	uint x = threadIdx.x;
	uint y = blockIdx.x;
	uint z = blockIdx.y;
	if ((x >= volumeSize.width) || (y >= volumeSize.height) || (z >= volumeSize.depth)) return;

	double _x = visualHullBox.x / volumeSize.width;
	double _y = visualHullBox.y / volumeSize.height;
	double _z = visualHullBox.z / volumeSize.depth;

	double a[2], b[2], c[2];

	double f[8][3];
	a[0] = visualHullBox.o[0] + _x * x;
	a[1] = a[0] + _x;
	b[0] = visualHullBox.o[1] + _y * y;
	b[1] = b[0] + _y;
	c[0] = visualHullBox.o[2] + _z * z;
	c[1] = c[0] + _z;
	int zz = 0;
	for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 2; j++) {
			for (int k = 0; k < 2; k++) {
				f[zz][0] = a[i];
				f[zz][1] = b[j];
				f[zz][2] = c[k];
				zz++;
			}
		}
	}

	unsigned short tag = 1;
	double minx, miny, maxx, maxy;
	int iminx, iminy, imaxx, imaxy;
	int ans;

	for (int i = 0; i < numberCamera; i++) {
		if (!tag) break;

		extrinsic_par[i].g_get_M_x_y(&minx,&maxx,&miny,&maxy,f);

		iminx = (int)minx;
		iminy = (int)miny;
		imaxx = (int)maxx + 1;
		imaxy = (int)maxy + 1;

		if (iminx < 0) iminx = 0;
		if (imaxx >= widthStep) imaxx = widthStep - 1;
		if (iminy < 0) iminy = 0;
		if (imaxy >= imageH) imaxy = imageH - 1;

		ans = integrogram[i * imageH * widthStep + iminy * widthStep + iminx] + integrogram[i * imageH * widthStep + imaxy * widthStep + imaxx] - integrogram[i * imageH * widthStep + imaxy * widthStep + iminx] - integrogram[i * imageH * widthStep + iminy * widthStep + imaxx];
		
		if (ans == 0) tag = 0;
	}

	if (tag) {
		visualHull[z * volumeSize.width * volumeSize.height + y * volumeSize.width + x] = 128;
	}else{
		visualHull[z * volumeSize.width * volumeSize.height + y * volumeSize.width + x] = 0;
	}

}
__global__ void
d_reconstruct_temperature(float* volumeData, uchar* visualHull,
	int numberCamera, int cameraIndex, int widthStep, int nChannels,
	char* imageData, uint *d_output, uint imageW, uint imageH,
	float density, float brightness,
	float transferOffset, float transferScale,
	int maxSteps, float tstep,
	float3 boxMin, float3 boxMax,
	cudaExtent volumeSize, unsigned long iteration)
{

	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;

	if ((x >= imageW) || (y >= imageH)) return;

	float4 color = IntToRgbaFloat(d_output[y * imageW + x]);
	float cur = color.y;


	float src = uchar(imageData[cameraIndex * imageH * widthStep
		+ (imageH - y - 1) * widthStep + x * nChannels + 1]);


	int threshold = 30;
	if (uchar(imageData[cameraIndex * imageH * widthStep
		+ (imageH - y - 1) * widthStep + x * nChannels + 0]) < threshold
		&& uchar(imageData[cameraIndex * imageH * widthStep
			+ (imageH - y - 1) * widthStep + x * nChannels + 1]) < threshold
		&& uchar(imageData[cameraIndex * imageH * widthStep
			+ (imageH - y - 1) * widthStep + x * nChannels + 2]) < threshold
		)
	{
		return;
	}
	float residual = src - cur;
	float lap = 0.f;
	//if (x < imageW - 1 && x > 0 && y < imageH - 1 && y > 0)
	//{

	//	float curUp = IntToRgbaFloat(d_output[(y - 1) * imageW + x]).y;
	//	float curDown = IntToRgbaFloat(d_output[(y + 1) * imageW + x]).y;
	//	float curLeft = IntToRgbaFloat(d_output[y * imageW + x - 1]).y;
	//	float curRight = IntToRgbaFloat(d_output[y * imageW + x + 1]).y;


	//	float srcUp = uchar(imageData[cameraIndex * imageH * widthStep
	//		+ (imageH - y) * widthStep + x * nChannels + 1]);
	//	float srcDown = uchar(imageData[cameraIndex * imageH * widthStep
	//		+ (imageH - y - 2) * widthStep + x * nChannels + 1]);
	//	float srcLeft = uchar(imageData[cameraIndex * imageH * widthStep
	//		+ (imageH - y - 1) * widthStep + (x - 1) * nChannels + 1]);
	//	float srcRight = uchar(imageData[cameraIndex * imageH * widthStep
	//		+ (imageH - y - 1) * widthStep + (x + 1) * nChannels + 1]);

	//	lap = (srcUp + srcDown + srcLeft + srcRight - 4 * src) - (curUp + curDown + curLeft + curRight - 4 * cur);
	//}

	if (fabsf(residual) < 10e-3)
	{
		return;
	}



	// 计算delat值完毕，后续进行强度调整
	// calculate eye ray in world space
	Ray eyeRay;
	eyeRay.o = c_nearFar[0] + x * (c_nearFar[3] - c_nearFar[0]) / imageW + y * (c_nearFar[1] - c_nearFar[0]) / imageH;
	eyeRay.d = c_nearFar[4] + x * (c_nearFar[7] - c_nearFar[4]) / imageW + y * (c_nearFar[5] - c_nearFar[4]) / imageH;
	eyeRay.d = normalize(eyeRay.d - eyeRay.o);

	// find intersection with box
	float tnear, tfar;
	int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);

	if (!hit) return;

	//if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

	// march along ray from front to back, accumulating color
	float t = tnear;
	float3 pos = eyeRay.o + eyeRay.d*tnear;
	float3 step = eyeRay.d*tstep;
	float3 normPos;
	float alpha = 1.f;
	float3 midPos = eyeRay.o + eyeRay.d * ((tnear + tfar) / 2.f);



	for (int i = 0; i<maxSteps; i++)
	{


		// read from 3D texture
		// remap position to [0, 1] coordinates
		// +0.5是由于cuda 3dtex的存储方式

		// 3d tex normalized coordinate
		normPos.x = (0.5f + ((pos.x) - boxMin.x) / (boxMax.x - boxMin.x) * (volumeSize.width - 1)) / volumeSize.width;
		normPos.y = (0.5f + ((pos.y) - boxMin.y) / (boxMax.y - boxMin.y) * (volumeSize.height - 1)) / volumeSize.height;
		normPos.z = (0.5f + ((pos.z) - boxMin.z) / (boxMax.z - boxMin.z) * (volumeSize.depth - 1)) / volumeSize.depth;

		float midDis = sqrt(dot(pos -  midPos, pos - midPos)) / (boxMax.x - boxMin.x);
		//printf("midDis = %f\n", midDis);

		// splatting
		//int scale = myRand() % 3 + 1;
		int scale = 1;
		int offset = 0;
		for (int xoffset = 0; xoffset < scale; ++xoffset)
		{
			for (int yoffset = 0; yoffset < scale; ++yoffset)
			{
				for (int zoffset = 0; zoffset < scale; ++zoffset)
				{
					d_adjust_temperature(volumeData, visualHull, normPos, volumeSize,
						xoffset - offset, yoffset - offset, zoffset - offset, midDis,
						residual, alpha, i, lap, src, numberCamera, iteration);
				}
			}
		}

		// "over" operator for front-to-back blending
		//alpha = alpha + density * ALPHA * (1.f - alpha);

		// "under" operator for front-to-back blending
		alpha = (1 - ALPHA) * alpha;


		t += tstep;

		if (t > tfar) break;

		pos += step;
	}

}


extern "C"
void reconstruct_temperature_kernel(dim3 gridSize, dim3 blockSize,
float* volumeData, uchar* visualHull,
int cameraIndex, uint *d_output, uint imageW, uint imageH,
float density, float brightness, float transferOffset, float transferScale,
int samplingMaxSteps, float samplingStep,
float3 boxMin, float3 boxMax,
cudaExtent volumeSize, unsigned long iteration)
{
	d_reconstruct_temperature << <gridSize, blockSize >> >(volumeData, visualHull,
		gcNumberCamera, cameraIndex, gSourceImageWidthStep, gSourceImageNChannels,
		d_source_image_data, d_output, imageW, imageH, density,
		brightness, transferOffset, transferScale,
		samplingMaxSteps, samplingStep,
		boxMin, boxMax, volumeSize, iteration);

	cudaThreadSynchronize();
}

extern "C"
void reconstruct_color_kernel(dim3 gridSize, dim3 blockSize, 
float* volumeData, uchar* visualHull,
int cameraIndex, uint *d_output, uint imageW, uint imageH,
float density, float brightness, float transferOffset, float transferScale, 
int samplingMaxSteps, float samplingStep,
float3 boxMin, float3 boxMax,
cudaExtent volumeSize, int channel)
{
	d_reconstruct_color<< <gridSize, blockSize >> >(volumeData, visualHull,
		gcNumberCamera, cameraIndex, gSourceImageWidthStep, gSourceImageNChannels,
		d_source_image_data, d_output, imageW, imageH, density,
		brightness, transferOffset, transferScale, 
		samplingMaxSteps, samplingStep,
		boxMin, boxMax, volumeSize, gChannel);

	cudaThreadSynchronize();
}

__global__ void d_smooth_visual_hull(VisualHullType* visualHull, cudaExtent volumeSize, int offset)
{
	uint x = threadIdx.x;
	uint y = blockIdx.x;
	uint z = blockIdx.y;
	//if ((x >= volumeSize.width) || (y >= volumeSize.height) || (z >= volumeSize.depth)) return;
	if (x >= volumeSize.width - offset
		|| y >= volumeSize.height - offset
		|| z >= volumeSize.depth - offset
		|| x < offset
		|| y < offset
		|| z < offset)
	{
		return;
	}

	unsigned long sum = 0;
	for (int xi = -offset; xi <= offset; ++xi)
	{
		for (int yi = -offset; yi <= offset; ++yi)
		{
			for (int zi = -offset; zi <= offset; ++zi)
			{
				sum += uint(visualHull[(z + zi) * volumeSize.width * volumeSize.height + (y + yi) * volumeSize.width + (x + xi)]);

			}
		}
	}


	visualHull[z * volumeSize.width * volumeSize.height + y * volumeSize.width + x] =
		VisualHullType(1.f * sum / (offset * 2 + 1) / (offset * 2 + 1) / (offset * 2 + 1));



}

extern "C"
void smooth_visual_hull(dim3 gridSize, dim3 blockSize,
	VisualHullType* visualHull,
	cudaExtent volumeSize, int offset)
{
	d_smooth_visual_hull<< <gridSize, blockSize >> >(visualHull,
		volumeSize, offset);

	cudaThreadSynchronize();
}

extern "C"
void generate_visual_hull_color_temperature_kernel_resize_octree(
	dim3 gridSize, dim3 blockSize,
	VisualHullType* visualHull,
	uint imageW, uint imageH,
	cudaExtent volumeSize, Octree* oct, uint oct_index, ResizeImage* gResizeImageData)
{
	d_generate_visual_hull_color_temperature_resize_octree << <gridSize, blockSize >> >(visualHull,
		gcNumberCamera, gSourceImageWidthStep, gSourceImageNChannels,
		gResizeImageData, imageW, imageH,
		volumeSize, visualHullBox, d_extrinsic_par, oct, oct_index);

	cudaThreadSynchronize();
}

extern "C"
void generate_visual_hull_color_temperature_kernel3(
	dim3 gridSize, dim3 blockSize,
	VisualHullType* visualHull,
	uint imageW, uint imageH,
	cudaExtent volumeSize, Octree* oct, Octree* oct_t,uint* quere, uint oct_index)
{
	d_generate_visual_hull_color_temperature3 << <gridSize, blockSize >> >(visualHull,
		gcNumberCamera, gSourceImageWidthStep, gSourceImageNChannels,
		d_source_imtegrogram, imageW, imageH,
		volumeSize, visualHullBox, d_extrinsic_par, oct, oct_t, quere, oct_index);

	cudaThreadSynchronize();
}

extern "C"
void generate_visual_hull_color_temperature_kernel2(
	dim3 gridSize, dim3 blockSize,
	VisualHullType* visualHull,
	uint imageW, uint imageH,
	cudaExtent volumeSize, Octree* oct, uint oct_index)
{
	d_generate_visual_hull_color_temperature2 << <gridSize, blockSize >> >(visualHull,
		gcNumberCamera, gSourceImageWidthStep, gSourceImageNChannels,
		d_source_imtegrogram, imageW, imageH,
		volumeSize, visualHullBox, d_extrinsic_par,oct,oct_index);

	cudaThreadSynchronize();
}

extern "C"
void generate_visual_hull_color_temperature_kernel(
dim3 gridSize, dim3 blockSize,
VisualHullType* visualHull,
uint imageW, uint imageH,
cudaExtent volumeSize)
{
	d_generate_visual_hull_color_temperature <<<gridSize, blockSize >>>(visualHull,
		gcNumberCamera, gSourceImageWidthStep, gSourceImageNChannels,
		d_source_imtegrogram, imageW, imageH,
		volumeSize, visualHullBox, d_extrinsic_par);

	cudaThreadSynchronize();
}
extern "C"
void reconstruct_color_temperature_kernel(dim3 gridSize, dim3 blockSize,
float* volumeData, uchar* visualHull,
int cameraIndex, uint *d_output, uint imageW, uint imageH,
float density, float brightness, float transferOffset, float transferScale,
int samplingMaxSteps, float samplingStep,
float3 boxMin, float3 boxMax,
cudaExtent volumeSize, int channel)
{
	d_reconstruct_color_temperature << <gridSize, blockSize >> >(volumeData, visualHull,
		gcNumberCamera, cameraIndex, gSourceImageWidthStep, gSourceImageNChannels,
		d_source_image_data, d_output, imageW, imageH, density,
		brightness, transferOffset, transferScale,
		samplingMaxSteps, samplingStep,
		boxMin, boxMax, volumeSize, gChannel);

	cudaThreadSynchronize();
}

__global__ void
d_compute_voxel_weight(short* pixelX, short* pixelY, int pixelSize, unsigned int* flameVoxelIndexArray, unsigned int flameVoxelSize, int camIdx, int camCount,
const uint imageW, const uint imageH,
float3 boxMin, float3 boxMax, float voxelLength,
int voxelStart, int voxelEnd,
uint3 voxelSize, int presetVoxelSize,
float* weight, unsigned int* indexFlameVoxel, unsigned int* weightSize)
{
	int indexPixel = blockDim.x * blockDim.y * blockDim.z * (blockIdx.z * gridDim.x * gridDim.y + blockIdx.y * gridDim.x + blockIdx.x)
		+ threadIdx.z * blockDim.x * blockDim.y + threadIdx.y * blockDim.x + threadIdx.x;

	if (indexPixel >= pixelSize)
	{
		return;
	}

	uint x = pixelX[indexPixel];
	uint y = pixelY[indexPixel];

	// calculate eye ray in world space
	Ray eyeRay;
	eyeRay.o = c_nearFar[0] + x * (c_nearFar[3] - c_nearFar[0]) / imageW + (imageH - y - 1) * (c_nearFar[1] - c_nearFar[0]) / imageH;
	eyeRay.d = c_nearFar[4] + x * (c_nearFar[7] - c_nearFar[4]) / imageW + (imageH - y - 1) * (c_nearFar[5] - c_nearFar[4]) / imageH;
	eyeRay.d = normalize(eyeRay.d - eyeRay.o);

	float tnear, tfar;
	int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
	if (!hit) return;

	uint3 indexVoxeluint3;
	float3 tempMinBox;
	float3 tempMaxBox;
	float tempNear, tempFar;

	for (int i = voxelStart; i < voxelEnd; ++i)
	{
		if (i >= flameVoxelSize)
		{
			break;
		}

		indexVoxeluint3 = calcPos(flameVoxelIndexArray[i], voxelSize);

		tempMinBox.x = indexVoxeluint3.x * voxelLength + boxMin.x;
		tempMinBox.y = indexVoxeluint3.y * voxelLength + boxMin.y;
		tempMinBox.z = indexVoxeluint3.z * voxelLength + boxMin.z;

		tempMaxBox.x = (indexVoxeluint3.x + 1) * voxelLength + boxMin.x;
		tempMaxBox.y = (indexVoxeluint3.y + 1) * voxelLength + boxMin.y;
		tempMaxBox.z = (indexVoxeluint3.z + 1) * voxelLength + boxMin.z;

		if (intersectBox(eyeRay, tempMinBox, tempMaxBox, &tempNear, &tempFar))
		{
			atomicAdd(weightSize + indexPixel, 1);
			indexFlameVoxel[indexPixel * presetVoxelSize + weightSize[indexPixel] - 1] = i;
			weight[indexPixel * presetVoxelSize + weightSize[indexPixel] - 1] = tempFar - tempNear;
			//weightSize[indexPixel] += 1;										
		}
	}
}
extern "C"
void compute_weight_kernel(short* pixelX, short* pixelY, const int pixelSize, unsigned int* deviceFlameVoxelIndex, unsigned int flameVoxelSize,
const int camIdx, const int camCount,
const uint imageW, const uint imageH,
const float3 boxMin, const float3 boxMax, const float voxelLength,
const int voxelStart, const int voxelEnd,
const uint3 voxelSize, const int presetVoxelSize,
float* weight, unsigned int* indexFlameVoxel, unsigned int* weightSize)
{
	int dimBlockX = 128;
	int dimGridX = (pixelSize / dimBlockX + 1 <= dimBlockX) ? (pixelSize / dimBlockX + 1) : dimBlockX;
	int dimGridY = (pixelSize <= dimBlockX * dimBlockX) ? 1 : (pixelSize / (dimBlockX * dimBlockX) + 1);
	dim3 dimGrid(dimGridX, dimGridY, 1);
	dim3 dimBlock(dimBlockX, 1, 1);
	//printf("blockx = %d, gridx = %d, gridy = %d\n", dimBlockX, dimGridX, dimGridY);

	if (dimGrid.x * dimGrid.y * dimGrid.z * dimBlock.x * dimBlock.y * dimBlock.z < pixelSize)
	{
		printf("not enough thread!\n");
		return;
	}



	d_compute_voxel_weight<< < dimGrid, dimBlock >> >(pixelX, pixelY, pixelSize, deviceFlameVoxelIndex, flameVoxelSize,
		camIdx, camCount,
		imageW, imageH,
		boxMin, boxMax, voxelLength,
		voxelStart, voxelEnd,
		voxelSize, presetVoxelSize,
		weight, indexFlameVoxel, weightSize);
}
#endif // #ifndef _VOLUMERENDER_KERNEL_CU_
